function result = opd_loop_wvl_a(input,fvindx)
    
    colvec = 'krgb'
    
    clear m_r_a m_i_a m_r m_i
    m_r = 1.48;
    m_i = 0.01*0.001*0.001;

    % particle log normal distribution
    % assume peak at 0.5 microns
    nx = 100;
    pd = makedist('Lognormal','mu',0.5e-4,'sigma',3.2);

    %r = linspace(0.1,1,nx)*1e-4;
    %r = [0.48e-4,0.49e-4,0.50e-4,0.51e-4,0.52e-4];
    r = [0.58e-4,0.59e-4,0.60e-4,0.61e-4,0.62e-4];
    nx = 5;
    %dens = pdf(pd,r);
    totdens = input/1.0; % dja scale 10x test 2025 01 06 
    dens = [1 1 1 1 1];

    %lognormal seems to be breaking
    %evenly distributing as a test
    dens = dens*totdens/(sum(dens))
    %dens = totdens;

    rho = 1.8; % g/cc (ish)

    % choose wavelength
    % hard coded kinetics grid
    wvls = [ 20.0, 40.0, 60.0, 80.0, 100.0, 120.0, 140.0, 160.0, 180.0, 200.0, ...
	       220.0, 240.0, 260.0, 280.0, 300.0, 303.8, 320.0, 340.0, 360.0, 368.1, 380.0, ...
	       400.0, 420.0, 440.0, 460.0, 480.0, 499.3, 500.0, 520.0, 540.0, 554.5, 560.0,...
	       580.0, 584.7, 600.0, 609.4, 620.0, 625.3, 629.7, 640.0, 660.0, 680.0, 700.0, ...
	       720.0, 740.0, 760.0, 764.6, 770.4, 780.0, 790.2, 800.0, 820.0, 834.2, 840.0, ...
	       860.0, 880.0, 900.0, 920.0, 923.1, 933.4, 937.8, 940.0, 949.7, 960.0, 972.5, 977.0, 980.0, ...
	       991.0, 1000.0, 1020.0, 1025.7, 1035.0, 1060.0, 1110.0, 1155.0, 1190.0, 1210.0, 1215.7, 1250.0, ...
	       1300.0, 1350.0, 1400.0, 1450.0, 1500.0, 1550.0, 1600.0, 1650.0, 1700.0, 1750.0, 1800.0, 1850.0, ...
	       1900.0, 1950.0, 2000.0, 2050.0, 2100.0, 2150.0, 2200.0, 2250.0, 2300.0, 2350.0, 2400.0, 2450.0, 2500.0, ...
	       2550.0, 2575.0, 2625.0, 2675.0, 2725.0, 2775.0, 2825.0, 2875.0, 2925.0, 2975.0, 3025.0, 3075.0, 3125.0, ...
	       3175.0, 3225.0, 3275.0, 3325.0, 3375.0, 3425.0, 3475.0, 3525.0, 3575.0, 3625.0, 3675.0, 3725.0, 3775.0, ...
	       3825.0, 3875.0, 3925.0, 3975.0, 4025.0, 4100.0, 4200.0, 4300.0, 4400.0, 4500.0, 4600.0, 4700.0, 4800.0, ...
	       4900.0, 5000.0, 5100.0, 5200.0, 5300.0, 5400.0, 5500.0, 5600.0, 5700.0, 5800.0, 5900.0, 6000.0, 6100.0, ...
	       6200.0, 6300.0, 6400.0, 6500.0, 6600.0, 6700.0, 6800.0, 6900.0, 7000.0, 7100.0, 7200.0, 7300.0, 7400.0, 7500.0, 7600.0, 7700.0, 7800.0, 7900.0, 8000.0]; %

    
    % hard coded a set of 4 as a test
    %wvls = [300.0,1000.0,6000.0,6100.0]; %
    wvls = wvls*1e-8;
  
    nwv = 175;
    %nwv=4;

    % Define the range of r
    min_r = r(1);  % Assuming r is sorted in ascending order
    max_r = r(end);

    % Find the indices corresponding to min_r and max_r
    idx_min = find(r >= min_r, 1);
    idx_max = find(r <= max_r, 1, 'last');

    %L = trapz(r(idx_min:idx_max), 33711083473.865845*3.141592*4.*r(idx_min:idx_max).^3/3.0) % part/cm2 x g/part = g/cm2
    L = totdens*1.8*3.141592*0.5e-4^3;

    tau = zeros(1,nwv);
    for iwv=1:nwv
	      x = 2*3.141592*r/wvls(iwv); %
              Qs = x*0;
              Qe = x*0;
              % number of terms in polynomial expansion
	      nmax = 100;
              n    = 1:nmax;
              a_n  = zeros(nx,nmax); %
	      b_n  = zeros(nx,nmax); %
	      m    = m_r + 1i*m_i; %
	      y    = m*x; %
	      
	      for j=1:nx % loop over Mie size parameter
		      [a, b] = mie_theory_fn(m,x(j),nmax);
	              a_n(j,:) = a; %
		      b_n(j,:) = b; %

		      % scattering and extinction efficiencies
		      Qs(j)    = sum((2*n+1).*(abs(a_n(j,:)).^2 + abs(b_n(j,:)).^2)); %
		      Qe(j)    = sum((2*n+1).*real(a_n(j,:) + b_n(j,:))); 
	      end
	      Qs = (2./x.^2).*Qs;
              Qe = (2./x.^2).*Qe;
	      beta = trapz(r(idx_min:idx_max), dens(idx_min:idx_max) .* Qe(idx_min:idx_max).* pi .* r(idx_min:idx_max).^2); %
	      ke = beta/(trapz(r(idx_min:idx_max), dens(idx_min:idx_max).*(rho).*(4/3.0).*3.141592.*r(idx_min:idx_max).^3));
	      tau(iwv) = ke*L;
    end

    a = 1-1e-4;
    mubar = 0.5;
    g = 0.85;
    rs = [0.49999896, 0.49998963, 0.4998963,  0.49976584, 0.49947142, 0.49880637, ...
	   0.49730681, 0.49391881, 0.48626963, 0.46900296, 0.43003111, 0.41703704, ...
	   0.40148148, 0.38592593, 0.36518519, 0.35573425, 0.35136986, 0.3439726, ...
	   0.30851507, 0.231332  ];
    rsrf = rs(fvindx);

    % delta scaling
    B = g^2;
    g_d = (g-B)/(1-B);
    tau_d = (1-B*a)*tau;
    a_d = ((1-B)/(1-B*a))*a;
    g = g_d;
    tau = tau_d;
    a = a_d;

    Gamma = 2*sqrt(1-a)*sqrt(1-a*g);
    r_inf = (sqrt(1-a*g) - sqrt(1-a))/(sqrt(1-a*g) + sqrt(1-a));

    denom = exp(Gamma.*tau) - r_inf^2.*exp(-Gamma.*tau);
    denom(100)

    r = tau;
    for i=1:nwv
	    r(i) = r_inf*(exp(Gamma*tau(i)) - exp(-Gamma*tau(i)))/denom(i);
    end
    
    r_inf
    tau(100)
    r(100)
    Gamma
    t_dir = exp(-tau/mubar);
    t_diff = (1-r_inf^2)./denom - t_dir;
    t = t_diff + t_dir;

    r_new = r
    for i=1:nwv	
	    r_new(i) = r(i) + rsrf*t(i)^2/(1-rsrf*r(i));
    end 
    r_new

    flux = [12800000.0, 34800000.0, 153000000.0, 153000000.0, 104000000.0, 53800000.0, 53800000.0, 737000000.0, 737000000.0, ...
                    553000000.0, 369000000.0, 369000000.0, 608000000.0, 436000000.0, 341000000.0, 6240000000.0, 346000000.0, 346000000.0, ...
                    84800000.0, 739000000.0, 84800000.0, 121000000.0, 157000000.0, 157000000.0, 303000000.0, 123000000.0, 8150000.0, 155000000.0, ...
                    203000000.0, 203000000.0, 799000000.0, 194000000.0, 194000000.0, 1580000000.0, 132000000.0, 450000000.0, 66300000.0, 3490000.0, ...
                    1500000000.0, 69800000.0, 88900000.0, 88900000.0, 469000000.0, 66700000.0, 66700000.0, 349000000.0, 200000000.0, 243000000.0, ...
                    349000000.0, 793000000.0, 561000000.0, 772000000.0, 38600000.0, 734000000.0, 1770000000.0, 1770000000.0, 1730000000.0, 1600000000.0,...
                    84300000.0, 84300000.0, 84300000.0, 1430000000.0, 84300000.0, 714000000.0, 35700000.0, 5960000000.0, 678000000.0, 54300000.0, 1030000000.0, 1460000000.0, 4380000000.0, 4460000000.0, 7560000000.0, 13500000000.0, 10800000000.0, 7860000000.0, 15600000000.0, 239000000000.0, 8080000000.0, 11900000000.0, 17900000000.0, 17100000000.0, 20000000000.0, 31100000000.0, 66000000000.0, 73300000000.0, 137000000000.0, 242000000000.0, 466000000000.0, 818000000000.0, 888000000000.0, 1680000000000.0, 2540000000000.0, 3490000000000.0, 5150000000000.0, 11900000000000.0, 18000000000000.0, 23400000000000.0, 29600000000000.0, 29500000000000.0, 26500000000000.0, 29900000000000.0, 36600000000000.0, 32300000000000.0, 32100000000000.0, 56600000000000.0, 97500000000000.0, 172000000000000.0, 142000000000000.0, 131000000000000.0, 173000000000000.0, 260000000000000.0, 417000000000000.0, 388000000000000.0, 397000000000000.0, 426000000000000.0, 502000000000000.0, 521000000000000.0, 596000000000000.0, 810000000000000.0, 838000000000000.0, 855000000000000.0, 869000000000000.0, 857000000000000.0, 892000000000000.0, 869000000000000.0, 947000000000000.0, 1080000000000000.0, 1030000000000000.0, 1040000000000000.0, 1000000000000000.0, 1040000000000000.0, 1060000000000000.0, 1310000000000000.0, 1700000000000000.0, 3680000000000000.0, 3840000000000000.0, 3570000000000000.0, 4100000000000000.0, 4610000000000000.0, 4740000000000000.0, 4790000000000000.0, 4870000000000000.0, 4780000000000000.0, 4870000000000000.0, 4880000000000000.0, 4760000000000000.0, 5020000000000000.0, 5020000000000000.0, 5080000000000000.0, 5060000000000000.0, 5200000000000000.0, 5360000000000000.0, 5320000000000000.0, 5320000000000000.0, 5300000000000000.0, 5220000000000000.0, 5240000000000000.0, 5240000000000000.0, 5090000000000000.0, 5120000000000000.0, 5230000000000000.0, 5220000000000000.0, 5110000000000000.0, 5140000000000000.0, 5010000000000000.0, 4930000000000000.0, 4910000000000000.0, 4810000000000000.0, 4810000000000000.0, 4750000000000000.0, 4640000000000000.0, 4650000000000000.0, 4610000000000000.0, 4550000000000000.0]    ;
    totflux = sum(flux)
    albedo = sum(r_new.*flux./totflux)
    
      result = albedo;
    end
  







